First principles theory of inelastic currents in a scanning tunneling microscope 



00 
On 
OS 

PL, 

l> 
(N 



K. Stokbro 1 , Ben Yu-Kuang Hu 1 , C. Thirstrup 2 and X. C. Xie 3 
1 Mikroelektronik Centret, Danmarks Tekniske Universitet, Bygning 3450, DK-2800 Lyngby, Denmark. 
2 Surface and Interface Laboratory, RIKEN, Saitama 351, Japan. 
3 Department of Physics, Oklahoma State University, Stillwater, OK 74078, U.S.A. 

(February 1, 2008) 

A first principles theory of inelastic tunneling between a model probe tip and an atom adsorbed on 
a surface is presented, extending the elastic tunneling theory of Tersoff and Hamann. The inelastic 
current is proportional to the change in the local density of states at the center of the tip due to the 
addition of the adsorbate. We use the theory to investigate the vibrational heating of an adsorbate 
below an STM tip. We calculate the desorption rate of H from Si(100)-H(2x 1) as function of the 
sample bias and tunnel current, and find excellent agreement with recent experimental data. 
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Atom manipulation using scanning tunneling micro- 
scopes (STMs) has been reported for a wide range of 
systems 0J^]. The manipulation of atoms has been at- 
tributed to various origins: direct forces between the tip 
and the surface |3|4|] , indirect interaction through the tip 
induced electric field |||| and inelastic scattering by elec- 
trons tunneling from the tip to the sur face |,|. The last 
mechanism, for example, is thought to account for the re- 
versible transfer of a Xe atom between a Ni(llO) surface 
and a W tip |P~|ll||, hydrogen desorption from hydrogen 
passivated silicon [Q, and dissociation of O2 molecules 
on a Ni(flO) surface ||. However, it is often not evident 
from the experimental data which mechanism dominates, 
and input from theoretical calculations is vital to estab- 
lish the microscopic mechanisms behind atomic manipu- 
lation on surfaces. 

In this paper we present a first principles method for 
calculating the inelastic current in a STM tunnel junc- 
tion due to inelastic scattering of the tunneling electrons 
with an adsorbate induced resonance. Our method is 
based on the Tersoff-Hamann approximation for the STM 
tunnel junction Jl2| and builds on previous work on in- 
elastic tunneling in STM The main result of this 
communication, Eq. (|^), is an expression for the inelas- 
tic current in terms of the partial local density of states 
(DOS) of the adsorbate wave function at the position 
of the tip. This expressions allows us to calculate the 
current and bias dependence of the desorption rate of H 
from the monohydrate Si(100) surface, by solving a Pauli 
master rate equation for transitions between the different 
vibrational levels of the adsorbate. In a recent letter |l6| ] 
we used Eq. (|^) to explain experimental desorption rates 
at negative sample bias due to inelastic hole tunneling. 
Here we present a derivation of Eq. ([)]), and use it to 
calculate desorption rates at positive sample bias due to 
inelastic electron tunneling. 

We divide the tunnel region into a tip and a sample 
part, and assume that the atom to be manipulated is 
adsorbed onto the sample. The model is shown schemat- 
ically in Fig. Jl| The total Hamiltonian is given by 
H = i?.5+a + -f/tip + Ht, where the terms on the right 



are the surface plus adsorbate, tip and tunneling parts, 
respectively. We first examine H^. Let \a) and |fc) be 
the one-electron eigenstates of the separated atom and 
the surface, respectively. When the atom adsorbs onto 
the surface, a coupling term h^a = Sfc tkac\a + h.c. is 
introduced (where Cfc and a are the field operators for 
states \k) and \a) respectively). The surface-l- adsorbate 
system can be diagonalized to yield new eigenstates 
The adsorbate also can vibrate on the surface, which im- 
plies a localized phonon term with a boson field opera- 
tor b. Furthermore, the adsorbate vibration couples to 
the energy-level of the electron in state \a). Therefore, 
H s +a = H e + i?ph + i? e _ p h, (the electron, phonon and 
coupling term respectively) where 



-Hph = b'b, 
He-ph = A(6 f + b)a)a, 



(la) 

(lb) 
(1c) 



The tip Hamiltonian is given by -f/tip = Y) p £pC P c p (where 
\p) are the tip eigenstates) and the tip to surface tunnel- 
ing by H T = J2 W Wj,Cp + h.c. = Y.kp t^pclcp + h - c - 
(here, k = k or a). The matrix elements are related via 

tup — JJ^wMHw 

The coupling between the tip and the surface is nor- 
mally weak, hence it is sufficient to work to lowest (i.e, 
second) order in the tunneling matrix clement. Let \p,n) 
be the state with an electron in state p on the tip and 
n phonons on the adsorbate, and be the coupled 

eigenstate of H s , arising from the interaction between 
electron state /j, and n' adsorbate phonons. The inelas- 
tic component of the current from tip to surface which 
changes the phonon occupation from n to n' is given by 
|7| 



HP 



\T P n !fm > \ 2 5(e p -s^- (n' - n)huj ), (2) 
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where T pn ^ pn i — (pn|-&r|/in') is the T-matrix element 
and / s ,t(e) = 1/{1 + exp[(e — £Fs,Ft)]} are Fermi distri- 
bution functions. The Fermi levels Ef s and ept of the 
sample and tip, respectively, are determined by the ap- 
plied sample bias, Vb, and related by Ef s = £_Ft — 

We define uncoupled electron-phonon eigenstates (i.e., 
setting A = 0) of the sample as ® |n). The coupled 
eigenstates can be expressed as a Born series H,f9i 



)H R - 



z=o 



I A 4 ) ® |n), 



(3) 



where e^. n « e p + n hajQ is the eigenenergy of n), and 
Gg. — (e — H e — /f p h + i0 + ) _1 is the Green function of 
the uncoupled electron-phonon system. Using Eq. (||) 
the T-matrix elements can be written as 



We now consider n' — n + N; i.e., where N phonons 
are absorbed (the emission case is similar) . The first term 
in the Born series Eq. (|J) which contributes is the Nth 
order term containing (&') . In this paper, we ignore the 
higher order terms in the series [EOl. This term gives 



where A(e) is the magnitude of the imaginary part of the 
self-energy of (a\g(e)\a). We now assume that the adsor- 
bate resonance is broad compared to the Hloq- Hence we 
(1) approximate the terms e^ k + (N— j)fkOo w E Pk , (2) as- 
sume A(e) « A(e„) = A, and (3) set ^ |(/i»| 2 S(e^~ 
e m ) w \(p,\a)\ 2 p s where p s = J^^ S(e-i a ) is the density of 
states at e a - These assumptions together with Eqs. (Q), 
(|) and (0) yield 



I(n -> n + N)-- 



{n + N)\2ire ^ N 



\(»\a)\ 2 » f t (e p )[l - f s (e,)]S(e p - Sfl - NTuj), (8) 

where K = irp s A 2 /A is a dimensionless quantity. 

The last step is to evaluate the matrix element t UhP — 
tkp using the Tersoff-Hamann approximation Jl2| , f2l| . 
^^Jp this approximation the tunneling is assumed to be 
through a single atom at the tip apex with an s-type 
wave-function. We define \t PkP - t kp \ 2 = / wp |f wp | 2 , 
where / Mfc = \{a\p)\ 2 / (x Pk + \(a\p,)\ 2 ) is essentially the 
fraction of electrons which tunnel from the tip into state 
Pk through the resonance. The evaluation of \t^ kP — tk p \ 2 
then is equivalent to previous work on clastic tunneling 
fl22fl and for a W tip we obtain 



TiZn+N = V(N + ny./n\ X N (o| M > £ t 



( a |ff( e /i,-ro + (N — j)huj Q )\a) (5) 
i=l,JV-i 

where g(e) = (e — _ff e + «0 + ) _1 is the sample+adsorbate 
electron Green function (without phonons). 

The Lippmann-Schwinger equation gives the one- 
electron eigenstates of the sample+adsorbate in terms 
of the eigenstates of the atom and the sample, \p k ) = 
\k) + g{e fj. k )\a)t* ka . This, together with Eq. (||) and the 
relationship between i w and t Kp yields 

T^Ln + N = V(n + N)!/n! \ N [t* pkp - tl 



J] ( a \g(e, k + (N-j)hu; )\a)+r^Z 

n at i 



j=0 s iV-l 



(6) 



The term T corresponds to processes where the electron 
tunnels from the tip into the sample state, k, hops back 
into the adsorbate, emits phonons and hops back out 
again. Generally, such terms are small, and moreover 
it may even be inaccurate to include this term, since the 
lifetime of state k probably will be dominated by scatter- 
ing with electron-hole excitations or bulk phonon states 
(which is excluded in our model). 

The next step is to relate the above equation to the 
projectors (a\p). We use the identity 



\(a\9(e)\a)\ 



2 _^IWI 2 <^ 



A(e) 



(7) 



I(n -> n + N) » C, 



(n + N)\ 
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Pn{<1 + i? w , e, E)de. 



(9) 



In this equation d is the tip-sample distance, k(e) = 
y/2m{<j)t + eVb — e)/h is the wave function inverse de- 
cay length, (fi t — 4.5 eV the work function of the W tip 
p3[ , E is the electric field between the tip and sample, 
i? w = 3bohr is the atomic radius of W, and the nor- 
malization constant C w = 0.007i?^ Amperes x bohr is 
obtained from a calculation of a model W tip p3|. The 
local weighted DOS, 

p n (r, e, E) = K n U\{a\p)\ 2n \M^ E)\ 2 S(e ~ e M ), 

is the DOS weighted by n powers of the dimensionless 
coupling K and the projection (a\fi) of the resonance 
wave function %fi a onto the eigenstates, ifi^, of the sample. 
For high biases ( | Vb | > 2 V) the effect of the electric field 
between the tip and the sample, E, must be taken into 
account when calculating ifi^ [ |22| . The inelastic current 
with energy transfer nhwo is weighted by (| (a\p) \ 2 K) n , 
thus only eigenstates with a significant overlap with the 
adsorbate resonance contribute to the inelastic current. 

We now use Eq. (Q) to model the STM induced des- 
orption of H from Si(100)-H(2x 1). It has previously 
been proposed [Q that inelastic scattering of the tun- 
neling electrons with the Si-H 6a* resonance is the domi- 
nant dcsorption mechanism for sample biases in the range 
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2 V < Vb < 4 V. In the following we make a first princi- 
ples calculation of the desorption rate due to this mecha- 
nism as a function of the tunnel current and sample bias 
and compare with experimental desorption rates. 

The electronic structure calculations are based on 
density-functional theory |^,^5| within the generalized- 
gradient approximation(GGA) pq] using 20 Ry plane- 
wave basis sets. We describe the Si(100)-H(2xl) surface 
by a 12 layer (2 x 1) slab and use 32 Appoints in the sur- 
face Brillouin zone. Ultra-soft pseudo-potentials [^7| are 
used for both H and Si. The geometry of the surface 
is obtained by relaxing the H atoms and the first six Si 
layers, and the resulting bond lengths and bond angles 
of the surface atoms are similar to other first principles 
calculations |2^,^|. From a frozen phonon calculation 
we obtain the Si-H stretch frequency tkoo = 0.26 eV, and 
subtracting the ground state energy of a free H atom we 
calculate the desorption barrier E^es = 3.36 eV. 

We will use the Si-H molecular 6a* eigenstate, |to6 CT *}, 
to represent the 6a* resonance. We expect \me a *) to 
have a large overlap with the 6cr* resonance, and smaller 
overlaps with the other resonances localized around Si-H, 
which are mainly the Si-H 4<r and 5cr resonances. The 
4<7 and 5cr resonances are well below the Fermi level, and 
we remove these components by projecting out the over- 
lap between \m 6a ^) and occupied slab eigenstates. In 
Fig. H we show the resulting 6a* partial DOS for dif- 
ferent Si-H separations, z. From the first and second 
moments we determine the average centers, £a(z), and 
widths, A(z), as shown in the inset of Fig. g. At the 
equilibrium separation, z — 0, the resonance is centered 
at e a — 4.6 eV, and has an average width A = 1 eV. For 
z > -0.1 A, the variation of e a with z is nearly linear and 
from a linear least squares fit in this region we calculate 
de a /dz — —3.12 ± 0.06 eV/A. The electron-phonon cou- 
pling term is given by A = y/h/ {2MuJo)de a /dz, where 
M is the adsorbate mass and from the above value of 
d£ a /dz we obtain A = 0.26 eV. 

To find p s we use the relation n a (e) w |(a|^)| 2 p s . Con- 
sidering only the eigenstates with the largest projections 
(a|/i), we find p s « 1.5 eV^ 1 . To estimate an average 
value of Xfj,, we will use the relation |i w | 2 = (|(£i|a)| 2 + 
Xn)\tap\ 2 - We first select the state, /3, with largest over- 
lap, (a\P), in the energy range of interest. If we assume 
\t ap \ an d x ii are constant in the energy range we obtain 
the equation x u = IWI'j^lM Wl 2 !*^! 2 which we uge 

to estimate x for each eigenstate in the energy range 2V- 
4V. We find the average value x « 0.01 (~ 0.1|(a|/3}| 2 ). 
With this value we obtain that Iq/I « 0.5 for biases in 
the range 2V-4V; thus 50% of the elastic electrons passes 
through the resonance. 

We now have all the parameters entering Eq. (^) for 
the inelastic current. From the inelastic currents we cal- 
culate the desorption rate, R, by solving the Pauli mas- 
ter equation for the transitions among the various levels 
of the oscillator JlC| ] , and vibrational relaxations due to 
phonon-phonon coupling is described by a current inde- 



pendent relaxation rate, 7 = 1 x 10 s s _1 0,|o). We 
assume that desorption occurs when the energy of the H 
atom exceeds the desorption energy Edcs = 3.36 eV, cor- 
responding to a truncated harmonic potential well with 
13 levels. Our objective is to obtain R as function of I 
and Vb- However, the inelastic current is a function of 
Vb and d, and we therefore calculate the elastic current, 
/(Vb, d), to establish a relation between /, Vb and d. For 
a given value of d and Vb we then calculate the inelastic 
currents, I(n — > n + N) and I(n + N — > n) for all bound 
vibrational states n = 1, 2, • ■ • 13. We include events with 
N = 1, 2, 3, and note that the N = 1 term gives the main 
contribution to the desorption rate in the current range 
relevant for desorption. 

In Fig. [| we compare the calculated current depen- 
dent desorption rate for sample biases of 2 V, 2.5 V, and 
3 V with experimental data obtained by us (open sym- 
bols) and Shen et al. {jj (filled symbols). The two exper- 
imental data sets have been obtained on highly doped n- 
type samples using similar experimental procedures. The 
dashed lines show the first principles calculation of the 
desorption rates. Both experimental and theoretical des- 
orption rates follow a power law R oc I a with a ~ 9 — 10 
for the theoretical curves and a ~ 10 — 13 for the experi- 
mental data. However, the theoretical curves are shifted 
towards lower currents. The discrepancy is due to the 
three main approximations of the theoretical model: an 
expansion based on A > A,wo, a harmonic approxima- 
tion for the Si-H bond potential, and neglect of excited 
state orbital relaxation when calculating the resonance 
wave function. In order to obtain theoretical desorption 
rates comparable to experiment we adjust the electron- 
phonon coupling. Solid lines in Fig. 3 shows results us- 
ing A = 0.20 eV, and we see that with this moderate 
change of the coupling constant we can obtain quantita- 
tive agreement with the experimental data in the voltage 
range 2-3 V. Above 3 V, the calculated desorption rates 
are too small to explain experimental desorption rates, 
and the measured desorption rates must be related to 
another mechanism, most likely direct excitation of the 
Si-H bond as suggested in Ref. Q . 

In conclusion, we have presented a novel first principles 
theory of inelastic scattering, and used it to calculate the 
voltage and current dependent variation of the hydrogen 
desorption rate from the Si(100)-H(2xl) surface at pos- 
itive bias conditions. We find that the desorption in the 
voltage range 2-3 V is consistent with vibrational heating 
of the H atom due to inelastic-scattering with the Si-H 
6a* resonance. 
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FIG. 1. Schematic figure of tunneling model (top) and asso- 
ciated density of states (bottom) . When the adsorbate (black 
dot) is separated from the surface, the eigenstates of the sur- 
face and adsorbate are |fc) and \a), respectively. Upon adsorp- 
tion, a coupling term with elements tka is introduced between 
|fc) and |a). In addition, there is a phonon term corresponding 
to the local oscillation of the adsorbate. The electronic part 
(i.e., when there are no phonons) can be diagonalized to yield 
eigenstates The center of resonance of the broadened 

level \a), e a , is a function of z the distance of the adsorbate 
from its equilibrium position. 
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FIG. 2. Projections of the Si-H molecular 6a* eigenstate 
onto the eigenstates, ft, of a Si(100)-H(2xl) slab for different 
Si-H separations(z = is the equilibrium separation). Only 
projections of states with e M > ef are included, and each par- 
tial DOS is normalized to 1. The inset shows the average 
centers(circle) and widths(error bar) as determined from first 
and second moments. The dashed line shows the results of a 
linear least squares fit to centers, e a , with z > — O.lA. 
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FIG. 3. Desorption rate, R, as function of current for 
Vt = 3.0 V(square), Vb = 2.5 V(circles) and 2.5 V(triangles). 
Open symbols show our experimental data, and filled symbols 
show data from Ref. . Dashed lines show theoretical calcu- 
lations with first principles value A = 0.26 eV for the electron 
phonon coupling, while solid lines show results using a fitted 
value A = 0.20 eV. 
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